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VESSEL SURFACE RECONSTRUCTION WITH A TUBULAR 
DEFORMABLE MODEL 

RELATED APPLICATIONS 

This application is based on, and claims benefit of, United States 
Provisional Application Serial Number 60/229,017 filed on August 30, 2000, 
and Untied States Provisional Application Serial Number 60/269,159 filed on 
February 1 5, 2001 . 

BACKGROUND OF THE INVENTION 

TECHNICAL FIELD OF THE INVENTION 
This invention generally relates to a data analysis method, apparatus 
and article of manufacture and more particularly to apparatus, article of 
manufacture and analysis method for analyzing arterial stenosis. 

DESCRIPTION OF RELATED ART 
Atherosclerosis is assessed according to the degree of stenosis of 
diseased vessels. For example, stenosis of the carotid artery of greater than 
70% represents a high risk of stroke and should be corrected by 
endarterectomy. The standard method for measuring the degree of stenosis 
is digital subtraction angiography (DSA) which provides high-contrast 
projection images of vessels from multiple viewpoints. However, DSA 
requires intra-arterial placement of a catheter that can have serious medical 
complications. Thus, alternatives to DSA have been proposed including 
magnetic resonance angiography (MRA), 3D ultrasound, and computed 
tomographic angiography (CTA), These methods are attractive because they 
are all less invasive than DSA and provide 3D images of the vessel lumen. 
However, the images from these methods are significantly more difficult to 
interpret due to lower image resolution and contrast. 

A variety of methods have been applied to the problem of vessel- 
surface reconstruction from 3D angiography. The simplest method for 
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reconstruction of vessel surfaces is iso-surface reconstruction. In this 
method, the location of the surface is determined strictly according to image 
intensity. A smooth surface is produced by placing the vertices at fractional 
image grid locations such that the interpolated image intensity is constant on 
the surface. However, choice of the threshold or iso-intensity level may be 
problematic, particularly for MRA, due to image inhomogeneities and other 
artifacts. Thus, the use of the iso-intensity surface is restricted to qualitative 
applications. 

Surface reconstruction of vessels has also been obtained from 
segmentation methods. Segmentation methods incorporate spatial and 
intensity criteria to objectively classify regions or voxels in an image. For 
example, a k-means clustering method automatically identifies voxels within 
vessels from MRA. This method compensates for the partial volume effect 
that diminishes the intensity of the smaller vessels thus causing errors for 
methods based on image intensity alone. Another method uses "fuzzy 
connectivity" for separating arteries and veins in MRA obtained with blood 
pool contrast agents. The arteries and veins are separated from the arteries 
based on a limited number of seed points inside and outside of the vessels. 
A similar watershed method has been successfully applied to MRA of the 
thoracic aorta. 

Once the image has been segmented, surfaces can be readily 
composed from segmentation boundaries although considerable smoothing 
must be applied and topological errors must be corrected. However, these 
segmentation methods incorporate only very limited a priori knowledge of the 
anatomy. Thus, while the shapes created by these methods directly reflect 
the information in the images, they tend to suffer unnecessarily from errors 
due to image artifacts or image noise. Several model-based approaches 
have been taken to address this problem. 

Several multi-scale methods have been applied to the problem of 
vessel surface reconstruction. In these methods, the width of objects in grey- 



scale images are automatically determined for any given point in the image 
by matching the objects with convolution kernels of various sizes. Thus, for 
points along a vessel, an average vessel diameter is obtained. Together with 
the vessel axis, the diameter measurements fully describe the vessel surface. 
5 However, all such multi-scale methods for vessel surface reconstruction 
assume that the vessel has a symmetric cross-sectional shape. Therefore, 
these methods have limited application for reconstruction of the vessel shape 
near bifurcations or in the presence of vascular disease such as 
atherorsclerosis. 

10 Another class of methods, the deformable models, allows for more 

general vessel shapes while still providing the useful smoothing effects. In 
the deformable models, surfaces deform according to image and smoothing 
"forces." The image forces pull the surface toward edges in the image while 
the smoothing forces resist bending of the surface and maintain even spacing 

15 of vertices of the surface mesh. Several deformable models have been 

applied to reconstruction of vessel surfaces. Two dimensional (2D) methods 
have been proposed in which the vessel boundary is reconstructed 
independently for each slice or plane. These provide adequate results for the 
cases tested but improvements can certainly be obtained by unified three 

20 dimensional (3D) reconstruction methods. 

Several general-purpose 3D deformable models have shown promise 
for surface reconstruction of large vessels. These methods have the 
potential to reconstruct non-tubular shapes such as vessel aneurisms and to 
provide an integrated reconstruction of bifurcations from minimal 

25 initializations. However, gross errors from these surface reconstructions 
have been found to occur in the region of bifurcations. 

Frangi et al. in "Model-based Quantization of 3-D Magnetic 
Angiographic Images" IEEE T. Med. Imaging, vol. 18, pp. 946-56 (1999), 
describe a 3D deformable model developed specifically for vessel surface 

30 reconstruction. The Frangi et al. model is based on a cylindrically 
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parameterized surface mesh. Since the cylindrical mesh is well-suited for 
vessel shape, a greater degree of smoothing can be obtained than from more 
generalized deformable model meshes. The model is initialized by a multi- 
scale method in which the axis is constructed while simultaneously estimating 

5 vessel diameter. Once the deformable model is initialized, the mesh deforms 
according to the conventional mechanical analogy. However, there are 
inherent limitations in this deformable model. One limitation is that the 
smoothing constraints used in this model inherently bias the surface 
reconstruction. For example, the "stretching" force on the surface vertices, 

10 which is meant to produce an even spacing of vertices, also tends to constrict 
the vessel. While a "bending" force is imposed to counteract that constriction 
effect, a complete cancellation cannot be obtained. Another inherent problem 
of this deformable model is the limited maximum density of vertices in the 
surface mesh. While this model can be readily initialized for a low density of 

15 vertices, surface self-intersection problems will occur if the density of vertices 
is increased in the axial direction. These shortcomings limit the accuracy of 
any vessel model, which may cause underestimation of vessel stenosis or, in 
the extreme, to completely overlooking the existence thereof. Further, 
measurement of vessel stenosis from 3D angiographic methods can be 

20 problematic due to limited image resolution and contrast. 

Thus, there is a need for improved, more accurate ways to display and 
analyze possibly diseased arteries. 

SUMMARY OF THE INVENTION 

It is a purpose of the present invention to improve reconstruction of 
25 vessel shape from 3D angiography; 

It is another purpose of the present invention to facilitate an objective 
evaluation of vessel shape; 

It is yet another purpose of the present invention to improve the 
precision of shape measurements from 3D angiography; 
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It is yet another purpose of the present invention to allow for more 
sophisticated characterization of vessel shape at stenoses; 

It is yet another purpose of the present invention to improve vessel 
surface reconstruction thereby facilitating accurate computational fluid 
dynamics (CFD) modeling of arterial blood flow; 

It is yet another purpose of the present invention to accurately 
estimate of fluid-mechanical conditions in a vessel that may be a factor in the 
progression of atherosclerosis. 

The present invention is an apparatus, article of manufacture, and 
method for modeling an elongated object located internal to a body. The 
elongated object may be a blood vessel such as the carotid artery or the 
renal artery. Magnetic resonance imaging data of the area of concern (e.g., 
the carotid artery) is collected. The magnetic resonance imaging data is 
analyzed, extracting gradient information. The gradient information may 
include the gradient of the magnitude gradient. Contemporaneously, a 
tubular coordinate system is interactively generated as an initial model of the 
artery. An axis and a reference circumferential direction are defined for the 
coordinate system with radial lines extending outward from the axis. 
Intersecting lines are merged. All vertices at radial and circumferential 
positions are initialized with the extracted gradient information. Then, the 
initialized model is deformed, subjecting initialized vertices to image and 
smoothing forces, thereby completing the surface model of the artery, 
effectively reconstructing the artery surface. The reconstructed artery 
surface may be displayed on a display. 

Thus, the tubular deformable model of the present invention is 
advantageous for quantification of vessel shape from 3D angiograph, 
particularly, where there is a mild rate of tapering of the lumen at a stenosis. 
In such subtle situations, the tubular deformable model of the present 
invention is likely to provide a more accurate measurement of the degree of 
stenosis than from manual interpretation. The preferred embodiment 
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deformable model provides enable more accurate realistic CFD modeling of 
blood flow patterns and better assessment of the risk of stroke, for example, 
posed by a given abnormality in the carotid lumen shape. 

In addition, the tubular deformable model of the present invention 
provides enhanced artery-vein separation in contrast-enhanced MR A, 
allowing application of the method of the present invention to a vein nearby to 
an artery of interest. Once the surface of the vein is reconstructed, it is 
nulled-out to leave an unobstructed view of the artery. 

The method is a deformable model that employs a tubular coordinate 
system. Vertex merging is incorporated into the coordinate system to 
maintain even vertex spacing and to avoid problems of self-intersection of the 
surface. 

BRIEF DESCRIPTION OF DRAWINGS 

The present invention will become more fully understood from the 
detailed description given hereinbelow and the accompanying drawings 
which are given by way of illustration only, and thus are not limitative of the 
present invention, and wherein: 

Figure 1 shows a schematic representation of a magnetic resonance 
analysis system for reconstructively modeling blood vessel surfaces with a 
tubular deformable mode! according to the preferred embodiment of the 
present invention. 

Figure 2 shows the reference circumferential orientation for a curved 

axis; 

Figure 3 shows an example of smoothing force deformation of surface 
mesh vertices about axis (a) analogous to mechanical deformation. 

Figure 4 shows an example of how points may be identified along the 
axis of the vessel. 

Figure 5 shows the reference circumferential orientation as defined for 
a curved axis tubular coordinate system. 
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Figure 6 shows a tubular coordinate system wherein the radial lines 
are warped in areas where the vessel axis is curved. 

Figure 7 shows meshes constructed from vertices at constant radial 
locations for a curved axis, i.e., on isodistance surfaces. 

Figure 8 shows surface reconstruction of phantom of carotid 
bifurcation reconstructed from an MR image of a glass phantom. A single 
slice of the image is shown in Panel A. Panels B and C show the surface 
reconstruction of overlapping segments. Panel D shows the bifurcation 
region by superimpostion of the two surfaces of B and C. 

Figure 9 is a example of surface reconstruction of carotid artery 
bifurcation from a normal subject. Panel A provides a cropped portion of the 
MRA in the sagital view. The segments are reconstructed separately (Panels 
B and C) and displayed as a superimposed image (Panel D) with shaped 
surfaces. 

Figure 10 shows surface reconstruction of performed on a carotid 
artery from a subject with stenosis of the carotid artery. Panel A provides a 
contrast-enhanced MRA. The segments are reconstructed separately in 
Panels B (external carotid artery) and C (internal carotid artery); the surface 
of the internal segment was manually indented at the grove between the 
internal and external carotid arteries (arrow in Panel C) where an error 
otherwise occurs due to discontinuity of the tubular shapes at that point. 
Panel D provides the superimposed image with shaped surfaces. 

Figure 1 1 shows surface reconstruction applied to MRA of patient with 
renal artery stenosis using the tubular deformable model. Panel A provides 
MRA. Reconstruction was performed independently on the aorta (Panel B) 
and the left renal artery (Panel C). The vessels are superimposed and 
visulaized as shaded surfaces in Panel D. 
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DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 

Turning now to the drawings and, more particularly, Figure 1 shows a 
schematic representation of a magnetic resonance analysis system 100 for 
reconstructively modeling blood vessel surfaces with a tubular deformable 
model according to the preferred embodiment of the present invention. It 
should be noted that, although the preferred embodiment is described herein 
with application to magnetic resonance imaging, it may be applied to 
reconstruction of images derived from other signal sources, e.g., acoustic 
data or acoustic resonance data, without departing from the spirit or scope of 
the invention. Thus, for example, the method of the present invention may be 
applied to ultrasound imaging. 

So, continuing, a pulsed radio frequency (RF) energy source 102 is 
directed to a patient 104 disposed in a main magnetic field generated by 
magnetic field generator 106. Receiver 108 receives magnetic energy 
emissions in response to the pulsed RF. Resonant magnetic energy signals 
received by the receiver 108 are passed to a analog-to-digital (A/D) converter 
110. Magnetic resonance digital data from the A/D converter 1 10 is passed 
to computer 1 12. Computer 112, converts the magnetic resonance digital 
data to viewable image data which is displayed on a suitable display screen, 
such as a cathode-ray tube 114, for example. 

Thus, according to the preferred embodiment of the present invention, 
image data collected from a patient 104 on the above system 100 is modeled, 
with allowances for curves in the vessel axis, variability in the vessel diameter 
and variability in cross-sectional shape. In the preferred deformable model, 
vessel shape is described by a radial function of axial and circumferential 
position, R(a, q>), Z 2 -»R within a tubular coordinate system. The tubular 
coordinate system, described in detail hereinbelow, resembles the cylindrical 
coordinate system but accommodates curved axes. The tubular coordinate 
system applied to a straight axis is identical to the well known cylindrical 
coordinate system. 
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Figure 2 shows an example of a flow diagram of the preferred 
embodiment method 120 for reconstructing vessel surfaces using a tubular 
deformable model. The preferred method 120 is bifurcated wherein gradient 
information is extracted in one leg 122 and, coincidentally, a tubular definition 
5 of the vessel is determined in the other leg 124. In step 126 the tubular 

information is combined with the gradient information, initializing the tubular 
vertices at all radial and circumferential positions. In step 128, the initialized 
tubular model is deformed, using additional gradient information, to result in 
the desired surface reconstruction. 

10 Thus, gradient information extraction begins in step 1222 as the image 

gradient is computed. Next, in step 1224 the magnitude of the gradient is 
computed. This gradient magnitude is provided to step 126 for initialization 
for the tubular model developed in step 124, such that all vertices fall at 
peaks in the gradient magnitude. Then, in step 1226 the gradient of the 

15 gradient magnitude is computed. This gradient of the gradient magnitude 
provides images and smoothing forces that are applied to the initialized 
shape in step 128, deforming the shape to complete reconstruction. 

Tubular shape definition begins in step 1242 wherein a user 
interactively defines the vessel being reconstructed. Preferably, this 

20 definition is done by taking slices through the vessel and the user identifying 
the vessel axis within the slice. Next, in step 1244 the vessel axis is 
interpolated using a b-spline method to connect the identified points. In step 
1246 the reference circumferential direction is defined. The reference 
circumferential direction is a function of axial position, a. Next, in step 1248, 

25 radial lines are defined for all axial and circumferential positions, a and q>, 

respectively. The radial lines extend outwards from the vessel axis. Tubular 
shape definition is completed in step 1250, when intersecting radial lines are 
merged. 

Accordingly, a discrete surface is described by a set of vertices that 
30 are evenly spaced in the axial and circumferential directions. Variation in 
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vertex position only occurs, however, in the radial direction. The radial lines 
are described as paths r a(p , constructed as provided hereinbelow. For a 
tubular coordinate system with a straight axis, the radial paths are straight 
lines projecting radially outwards from the axis. 
5 So, in step 1222, VI is obtained by convolution of the image with the 

gradient of the normalized spherical gaussian function (space constant, o = 1 
voxel unit). The radial function R(a, <p) is the value of the radial parameter r 0 
at the maxima of the gradient magnitude. Next in step 1224, the magnitude 
of the gradient is computed such that: 

10 IIV/G^))! > HV/(/^(r))|| VrE (0,r limit ) (1) 



where r,i mH is the value of the radial parameter at the furthest extent of the 
radial line. A further provision is that a decrease in image intensity occurs at 
that point proceeding outwards from the vessel center: 



V/(p^(r 0 ))-v a /r 0 ) > 0 (2) 
15 where v (r) is a unit vector in the direction of the radial line: 



\\p' a9 ir)W 



(3) 



The radial function deforms so as to minimize discontinuities in the 
radial position along the surface while maximizing the proximity of the 
vertices to edges in the image. The result is an equilibration process 
20 analogous to that of a mechanical system acted on by internal elastic forces 
and external forces. 

In step 1226, the gradient of the previously computed gradient 
magnitude is found for use in determining the instantaneous deformation of 
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vertices. The instantaneous deformation of a vertex at any axial- 
circumferential location is given by: 



Ar = K l (K 2 X^^-^ + v^C^-VUV/C^Ci?))!!) (4) 

neighbor 

where R nei ghbor is tne radial location of the neighbors of a given vertex within 
the surface mesh. The deformations are applied simultaneously to all 
vertices in step 126 and repeated until an equilibrium condition is obtained in 
step 128. An example of suitable source code for carrying out the present 
method is attached as an Appendix to United States Provisional Application 
Serial Number 60/269,159 (filed February 15, 2001); this source code, which 
is hereby incorporated by reference, is designed to run as part of IRIS 
Explorer platform on a Silicon Graphics Unix workstation. 

Figure 3 shows an example of smoothing force deformation of surface 
mesh vertices 130, 132, 134 about axis (a) 136, analogous to mechanical 
deformation. All surface vertices 130, 132, 134 deform along radial lines, r a(p , 
130. The negative gradient magnitude image provides a gravitational force 
(curved line 139) analog on the vertex 132, while the neighboring vertices 
pull the vertex with spring-like forces. The spring-like smoothing forces are 
proportional to the difference in the radial locations between the vertices 130, 
132, 134. 

Additional smoothing forces may be included in the deformation 
process. For example, forces may be applied through user interaction to 
correct for gross errors in the surface reconstruction. An "anchor" manually 
placed at location (a manual , q> manua|) r manua ,) will attract the surface towards itself. 
The modified deformation is given by: 



^manual = & + , (R - ) 

where S{a manual i q> nmtua ) = 1 for (a, <p) = 



fa manual' *Pmanuat) 



(5) 
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= 0 for (a, q>) * (a manual , <p manu J. 



As noted hereinabove, the preferred embodiment deformable model of 
the vessels is based on a tubular coordinate system that allows for curvature 
of the cylindrical axis. The tubular coordinate system differs from the 
cylindrical coordinate system in several ways. Accordingly, the tubular 
coordinate system is constructed from a vessel axis. The axis is formed from 
a sequence of points along the center of the vessel identified by the user in 
step 1242. A sufficient number of points must be chosen to represent all the 
curves in the vessel. Points along the center of the vessel can be easily 
identified by the user by displaying slices of the image normal to the axis of 
the vessel side-by-side with a 3D surface display of the vessels to provide 
anatomical reference. 

Figure 4 shows an example of how points may be identified along the 
axis of the vessel 1242. A cropped portion of the carotid MR A is shown by 
itself 140 and with a 3D isosurface 142 to provide anatomical context that 
helps to follow a given vessel from one axial slice to another. In certain 
cases where the vessei follows a torturous path, points along the axis may be 
chosen from the maximum intensity projection (MIP). Identification of points 
on a MIP provides a 3D coordinate since each point in the MIP is associated 
with a single point in the 3D image in the projection process. In step 1244, 
the axis is smoothed with a cubic b-spline and interpolated to have a spacing 
of 0.5 mm between points. 

Figure 5 shows the reference circumferential orientation as defined in 
step 1244 for a curved axis tubular coordinate system 150, as represented by 
raised ridge 152. The reference angle on which the circumferential angle is 
based is dependent on the curvature of the vessel axis. The circumferential 
orientation cannot be defined with respect to an absolute orientation due to 
the curvature of the axis. Rather, the circumferential component of the 
tubular coordinate system must be defined relative to the local curvature of 
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the axis. A circumferential coordinate is defined such that a minimum of 
angular change occurs between adjacent points on the axis. Taking A(a) as 
the path of the vessel axis, ju(a) as the unit vector along the vessel axis at 
the appoint, andv 0 (c?) as the reference radial orientation at the a th axial point: 

^(a)^o(0)- (7) 
The reference orientation at the adjacent point is then: 



v ( a ii)- -/< fl+1 X(V g )^ g+1 )) 

\\*M -Ka+mM-M(a+l))\\' ( ) 



A normal radial vector is obtained by the cross product of the 
reference radial vector with the axis vector: 

v v(«)= v » x ^). (9) 

All initial radial vectors are constructed from a linear combination of the 
reference and normal radial directions. Sixteen (16) circumferential 
directions of each axial location are used, evenly circumferentially spaced 
about the axis. Thus: 

v(a,4,0) = v w (a). (10) 
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A second feature of the preferred tubular coordinate system is that the 
radii do not emanate in straight lines from the axis at all points. As can be 
seen from the tubular coordinate system in the example of Figure 6, the radial 
lines are warped in areas where the vessel axis 160 is curved. Radial lines 
5 from an axis 160 will intersect at the inside of axis curves, i.e., in area 162. 
This problem is solved in step 1250 by merging radial lines in the tubular 
coordinate system. Radial lines are shown for the in-plane radial directions 
along a curved portion of the axis of a tubular coordinate system. The radial 
coordinate at each point in the tubular coordinate system is the distance to 
10 the nearest point on the axis. This warping prevents radial lines from adjacent 
axial locations from intersecting one another. 

The radial lines are defined prior to the deformation process, i.e., in 
step 1248. The radial lines extend outwards at a constant step size 
(Step_size) of 0.1 of the in-plane pixel resolution: 

1 5 p a<p (r) = A(a) + v(a, (p,0) xrx Step_size (11) 

and in step 1250, radial lines are concatenated when intersections would 
occur. In principle, the radial lines from an axis will not intersect one another 
exactly due to the discrete locations of the radial lines and to out-of-piane 
components of the radial direction vectors. Thus, practical definitions of line 
20 intersections must be adopted. The radial lines are considered to intersect 
one another in one of two ways: 

(1) A radial line intersects with an oblique or anti-parallel radial line if 
it reaches a point where it is closer to the origin of the second radial line than 
it is to its own origin. The intersection point of the radial line is designated, 

\\p av {r)-A{a)\\ < \\p a ^r)-A(a')\\ 
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V r < r franc and |«-a'|>l (12) 

(2) A radial line intersects with a nearly parallel line when it is at its 
closest approach to the second line. The point of intersection on the second 
line is designated r int 

where: 

WPa^^)- p a , 9 {r mt )\\ < \\p a ,(r)-p a , p (r')\\ V {rj') and \ a-a'\ =1. (13) 

The second condition applies specifically to the intersection between 
two radial lines from adjacent points on the axis while the first condition 
applies to all other intersections. These criteria are only applied to radial 
lines at the same circumferential positions. The truncation of any given radial 
line is the minimum of r< runc obtained from equations 12 and 13. 

The radial lines are then modified by merging so that they extend 
outwards to the maximal radial position. 
For r = r ]im „-l to r = 1 

For all (a,<p) 

if r tr l mc( a ^) = r 

for r'-r to r ]imit 

arrest the axial location of the non-concatenated radial line that is nearest 
in the axial direction such that: 
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unMn e arest^) = r imlt . (14) 



Figure 7 shows meshes 170, 172, 174 constructed from vertices at 
constant radial locations for a curved axis, i.e., on isodistance surfaces. 
Surfaces are constructed whose vertices are at equal distances from the 
5 given axis. The meshing is done according to the tubular coordinate system. 
For small radii mesh 170, vertices exist for each circumferential and axial 
position. At larger radii meshes 172, 174, the vertices merge with one 
another in the axial direction to eliminate criss-crossing of radial lines. 

The cylindrical deformation process could be applied directly to the 

10 modified cylindrical coordinate system. However, bunching of vertices will 
occur where radial lines merge. Where such bunching occurs, the effective 
elasticity of the surface becomes significantly greater which reduces the 
surface smoothness. This effect is removed by merging vertices in the 
deformation process to match the merging of the radial lines in step 1248. 

15 This precludes two vertices from co-existing on the same radial line. The 

procedure for merging vertices is as follows. The vertices initially exist for all 
circumferential and axial locations. After a primary deformation phase, 
vertices are merged with their axial neighbors if either of the vertices is at a 
point where the radial lines of the two vertices are merged: 

20 Merge (a,<p) and {a adjacent ,y) if p a9 (R) = p amtac ^(R) . 

The radial location of the resulting vertex becomes the average of the 
locations of the vertices from which it was formed. 

R (<*n,er ged , <P) = ~ ~ ( 1 5) 
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where N vertices are merged. 

Henceforth this vertex deforms as a single vertex. The neighbors of 
this resultant vertex are those of the merged vertices. However, the 
circumferential neighbors for the resultant vertex will accumulate relative to 
the axial neighbors. Therefore, to maintain a balance between the axial and 
circumferential component of the elastic force, the total circumferential elastic 
force must be rescaled. With regard to the equation of motion 4 the elastic 
component R - R neighbon the difference of radial location between a vertex and 
its neighbors, must be modified so that the circumferential and axial radial 
differences are added separately. So, 

R-R^ =i(x-Zj)+jft(R-R k ) (16) 

where K is the number of circumferential neighbors of a merged vertex. 

Figure 8 shows surface reconstruction of phantom 1 80 of carotid 
bifurcation. The carotid artery bifurcation is reconstructed from an MR image 
of a glass phantom. A single slice 180 of the image is shown (Panel A). The 
surface was reconstructed with the tubular deformable model of overlapping 
segments 182, 184 (Panels B and C). The superimposition of the two 
surfaces accurately represents the bifurcation region 1 86 (Panel D). The 
tubular deformable model was applied to an MR image of a glass phantom of 
the carotid artery bifurcation with stenosis of the internal carotid artery (CA) 
(3D Spoiled Gradient Recalled Echo (SPGR), Voxel dimensions 
0.31x0.31x1. 0mm, Phantom filled with 60:40 ratio water to glycerine, No fluid 
flow). The glass phantom was constructed according to typical dimensions 
(common CA (CCA)=8mm, internal CA (ICA)=7mm, external CA (ECA)=6mm) 
with a moderate stenosis at the origin of the ICA. The CCA/ICA segment and 
CCA/ECA segments of the phantom were independently reconstructed from 
the MR image. The same value of the elasticity parameter, K 2 = 7x1 0 3 cm" 1 is 



- 17- 



Docket 71527 



10 



used for both the phantom and data from human subjects. The space 
constant of tne gradient magnitude operator, vi (a, <p,r) , is 1.0 in-plane voxel 
units for all tests of the algorithm. Image contrast is similar for all images. 
The contrast is in the range of 200-300 gray-scale units. The reconstruction 
of the phantom is carried out in a total of one hundred deformation iterations. 
Fifty iterations are carried out prior to merging of vertices and fifty after 
merging of vertices. The iteration step size, K, was 0.05 cm. 

An average vessel radii was obtained at each axial position along the 
vessel segments. To calculate the average radius, a cross-sectional area 
was calculated from the vertex positions at that radial position: 



where T ( - is the triangular area between the axial point and two 
circumferential ly adjacent points and P } is the polygonal area at the j th axial 
location. The average radius is then given by the formula for the area of a 
circle multiplied by a constant factor to compensate for the difference in area 
between a triangle and a corresponding arc of a circle 



The average radii of a segment was calculated from 8-10 axial points. 
The average radii of all segments differed from the known radii by less than 
20% of the average voxel dimension. Radii at the stenosis were compared 
with radii measured manually from computed tomography image (CT) of the 
phantom (15cm FOV, reconstructed to 0.3 mm in-plane resolution, oriented 
normal to axis of stenosis). The difference between the radii measured by 
the deformable model and from CT was 0.2 mm. 
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Figure 9 is an example of surface reconstruction of carotid artery 1 90 
bifurcation from a normal subject. A cropped portion of the MRA is visualized 
as a MIP in the sagital view 192 (Panel A). The CCA/ECA segment 194 
(Panel B) and the CCA/ICA segment 196 (Pane! C) are reconstructed 
5 separately and superimposed as a shaded surface display 198 (Panel D). 
Surface reconstruction of the carotid artery 1 90 was carried out on contrast- 
enhanced MRA images. MR images were acquired with contrast injection of 
Gd-DTPA (fast-gradient recalled echo sequence, GE Corp.) of one carotid 
artery from each of 4 normal subjects. As with reconstruction of the phantom 

10 180, the reconstruction was carried out with fifty deformations prior to 

merging of vertices and fifty after merging. The iteration step size was also 
the same as for the phantom reconstruction (i.e., K-, = 0.05 cm). 

As with the example of Figure 9, realistic surfaces were obtained for all 
4 images that were consistent with both the MIP's and the source images. 

15 Only a limited portion of ECA and its sub-tree was reconstructed since the 
vessel rapidly bifurcates which increases the complexity of identifying the 
axis. However, the portion of the ECA reconstructed is adequate for CFD 
modeling of fiow in the carotid bifurcation. 

Reconstruction was performed on each of two carotid arteries from 

20 each of two subjects with stenosis of the carotid artery of which, Figure 1 0 
shows surface reconstruction of one, reconstructed with the tubular 
deformable model according to the preferred embodiment of the present 
invention. A carotid artery bifurcation from contrast-enhanced MRA, shown in 
MIP 200 (Panel A), is reconstructed as one segment continued onto the ECA 

25 202 (Panel B) and one onto the ICA 204 (Panel C). The two surfaces are 
superimposed in 206 (Panel D). The surface of the internal segment was 
manually indented at the groove between the internal and external carotid 
arteries (arrow 208 in Panel C) where an error otherwise occurs due to the 
discontinuity of the tubular shape at that point. In another example (not 

30 shown), the reconstructed surface was found to be smooth and entirely 
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consistent with the MIP and source images. For the image shown in Figure 
10, the surface reconstruction has one noticeable error at the origin of the 
ICA. At this location, an indentation in the lumen was smoothed-over in the 
reconstruction. This error occurred due to small size of the indentation and 
5 due to proximity to the insertion of the ECA where there is a discontinuity in 
the radial location of the surface. Such an error can be manually corrected 
by inserting a point to provide an attractive force towards the correct location 
as described hereinabove. The elastic coefficient of the manual attractive 
force K 3 was 7x1 0 4 cm" 1 (in equation 4). When a manual attractive force is 

10 present the iteration step size was reduced to K, = 0.025 cm to avoid 

instability in the simulation of the mechanical deformation. The total number 
of deformation iterations is then increased to a total of one hundred. 

The degree of stenosis apparent in the surface reconstruction is less 
than in the MIP. However, the discrepancy is not due purely to smoothing 

15 effects of the deformable model since a similar degree of stenosis is 

observed when the surface is constructed with no smoothing. The error may 
be attributed to the tendency of the MIP to over-estimate the degree of 
stenosis. Only one artery from each subject was reconstructed since, in the 
above example of Figure 10, the contra-lateral vessel was occluded in one 

20 case and suffered from severe image artifacts in the other. 

Figure 1 1 shows surface reconstruction applied to an MRA 210 (Panel 
A) of patient with renal artery stenosis using the tubular deformable model 
1 10 according to the preferred embodiment of the present invention. 
Reconstruction was performed independently on aorta 212 (Panel B) and on 

25 left renal artery 214 (Panel C). These vessel reconstructions are 

superimposed and visualized as shaded surface 216 (Panel D). For this 
example, the abdominal aorta and renal artery were reconstructed in two 
subjects with renal artery stenosis from contrast-enhanced MRA (fast 
gradient recalled echo sequence, GE Corp.). The surface reconstruction 

30 parameters were the same as those used for the carotid arteries. Surfaces of 
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both the renal arteries and the aortas are smooth and are consistent with MIP 
and source images. 

Thus, the tubular deformable model of the present invention provides a 
more attractive method than other state of the art methods for blood-vessel 
5 surface reconstruction from 3D angiography. While the Frangi et al. model 
for vessel surface reconstruction, for example, also employs a cylindrical 
surface mesh, the Frangi et al. model is quite different from the model of the 
present invention in important respects. One important advantage of the 
method of the present invention is found in initialization of the deformable 

10 model in step 126. The multi-scale method used by Frangi et al. (and also 
described by Aylward et al. in "Intensity Ridge and Widths for Tubular Object 
Segmentation and Description," Proceedings of IEEE Workshop, 
Mathematical Methods in Biomedical Image Analysis, June, 1996) requires a 
consistent cross-sectional image intensity profile along blood vessels. 

15 However, for a contrast-enhanced MR A, flow and timing artifacts can 

significantly distort the cross-sectional image intensity profile. For example, 
the jugular vein, which is immediately adjacent to the carotid artery will 
enhance to varying degrees, depending on the rates of venous return in the 
given individual and to the exact timing of the image acquisition. The 

20 presence of the jugular vein next to carotid artery affects the localization of 
the axis of an object using Frangi et al. methods. By contrast, initialization of 
the tubular deformable model is a one dimensional (1 D) edge detection 
method which assumes only that the image intensity decrease is greatest at 
the boundary of the lumen. 

25 Another advantage of the method of the present invention is that 

meshes are constructed with high vertex density. This allows for 
representation of smaller-scale detail of the lumen shape such as focal 
stensoses or protrusions into or out from the lumen. High vertex density 
cannot be obtained for simple cylindrical meshes where vertices are simply 

30 placed along normals from the axis as in Frangi et al. For such a cylindrical 
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mesh, the radial lines from nearby points along the axis will converge at the 
inside of curves of the cylindrical axis and intersect with one another. 
Advantageously, this problem is avoided in the tubular coordinate system of 
the present invention by a systematic merging of radial lines. 
5 Further, the preferred tubular coordinate system allows for the 

description of the vessel shape in terms of a single radial parametric function. 
This has a subtle advantage over the rectilinear parametric functions used in 
the model of Frangi et al. wherein smoothing of rectilinear parametric shape 
functions must be biased towards either under- or over-estimation of vessel 
10 diameter. By contrast, no such bias is present in the tubular deformable 
model of the present invention. While the rectilinear parametric shape 
functions of Frangi et al. allow for reconstruction of more generalized shapes, 
only a modest advantage is realized for the reconstruction of blood vessel 
surfaces. 

15 Thus, the tubular deformable model of the present invention is 

advantageous for quantification of vessel shape from 3D angiograph, 
particularly, where there is a mild rate of tapering of the lumen at a stenosis. 
In such subtle situations, the tubular deformable model of the present 
invention is likely to provide a more accurate measurement of the degree of 

20 stenosis than from manual interpretation. The preferred embodiment 

deformable model provides more accurate realistic CFD modeling of blood 
flow patterns and better assessment of the risk of stroke, for example, posed 
by a given abnormality in the carotid lumen shape. Once the vessel surface 
has been reconstructed, more sophisticated characterization maybe used to 

25 characterize vessel shape at stenosis. So, for example, the degree of 
stenosis may be easily derived by comparing the typical diameter of the 
vessel with the diameter at the stenosis and, in response, automatically 
suggesting a course of treatment, e.g., surgery and/or medication. 

In addition, the tubular deformable model of the present invention 

30 provides enhanced artery-vein separation in contrast-enhanced MR A, 
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allowing application of the method of the present invention to a vein nearby to 
an artery of interest. Once the surface of the vein is reconstructed, it is 
nulled-out to leave an unobstructed view of the artery. 

The method is a deformable model that employs a tubular coordinate 
system. Vertex merging is incorporated into the coordinate system to 
maintain even vertex spacing and to avoid problems of self-intersection of the 
surface. The deformable model was evaluated on clinical magnetic 
resonance (MR) images of the carotid (n=6) and renal (n=2) arteries and on 
an MR image of a vascular phantom. Only one gross error occurred for all 
clinical images. All reconstructed surfaces had a realistic, smooth 
appearance. For all segments of the phantom, vessel radii from the surface 
reconstruction had an error of less than 0.2 of the average voxel dimension. 

The invention being thus described, it will be obvious that the same 
may be varied in many ways. Such variations are not to be regarded as a 
departure from the spirit and scope of the invention, and all such 
modifications as would be obvious to one skilled in the art are intended to be 
included within the scope of the following claims. 
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